/*
 * OREKIT-X
 * Copyright 2002-2008 CS Communication & Systemes
 * 
 * Parts of this software package have been licensed to CS
 * Communication & Systemes (CS) under one or more contributor license
 * agreements.  See the NOTICE file distributed with this work for
 * additional information.
 *  
 * This is an experimental copy of OREKIT from www.orekit.org.
 * Please use the original OREKIT from orekit.org for normal work
 * unrelated to this research project.
 * 
 * Licensed under the Apache License, Version 2.0 (the "License"); you
 * may not use this file except in compliance with the License.  You
 * may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
 * implied.  See the License for the specific language governing
 * permissions and limitations under the License.
 */
package ore.forces.radiation;

import ore.CelestialBody;
import ore.ForceModel;
import ore.Frame;
import ore.PositionVelocity;
import ore.SpacecraftState;
import ore.errors.OrekitException;
import ore.propagation.events.AbstractDetector;
import ore.propagation.events.EventDetector;
import ore.propagation.numerical.TimeDerivativesEquations;
import ore.time.AbsoluteDate;

import org.apache.commons.math.geometry.Vector3D;

/**
 * Solar radiation pressure force model for a spacecraft in orbit
 * about a body in its reference frame.
 *
 * @author Fabien Maussion
 * @author &Eacute;douard Delente
 * @author V&eacute;ronique Pommier-Maurussane
 */
public class SolarRadiationPressure
    extends Object
    implements ForceModel
{

    private static final String LOW_TRAJECTORY_MESSAGE =
        "trajectory inside the Brillouin sphere (r = {0})";


    private static final double SUN_RADIUS = 6.95e8;

    /**
     * Reference distance (m).
     */
    private final double dRef;
    /** 
     * Reference radiation pressure at dRef (N/m<sup>2</sup>).
     */
    private final double pRef;
    /** 
     * Sun model.
     */
    private final CelestialBody sun;
    /** 
     * Earth model.
     */
    private final double equatorialRadius;
    /** 
     * Spacecraft.
     */
    private final RadiationSensitive spacecraft;


    /** Default reference values.
     * <ul>
     *   <li>d<sub>ref</sub> = 149597870000.0 m</li>
     *   <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m<sup>2</sup></li>
     * </ul>
     * @param sun Sun model
     * @param equatorialRadius Spherical shape model for umbra/penumbra computation
     * @param spacecraft  Physical and geometric information
     */
    public SolarRadiationPressure(CelestialBody sun, double equatorialRadius,
                                  RadiationSensitive spacecraft)
    {
        this(149597870000.0, 4.56e-6, sun, equatorialRadius, spacecraft);
    }

    /** 
     * @param dRef Reference distance for the radiation pressure (m)
     * @param pRef Reference radiation pressure at dRef (N/m<sup>2</sup>)
     * @param sun Sun model
     * @param equatorialRadius Spherical shape model for umbra/penumbra computation
     * @param spacecraft Object physical and geometrical information
     */
    public SolarRadiationPressure(double dRef, double pRef,
                                  CelestialBody sun,
                                  double equatorialRadius,
                                  RadiationSensitive spacecraft) {
        this.dRef  = dRef;
        this.pRef  = pRef;
        this.sun   = sun;
        this.equatorialRadius = equatorialRadius;
        this.spacecraft = spacecraft;
    }

    /** {@inheritDoc} */
    public void addContribution(SpacecraftState s, TimeDerivativesEquations adder)
        throws OrekitException
    {
        /*
         * Raw radiation pressure
         */
        Vector3D sunSatVector = getSatSunVector(s).negate();
        final double dRatio = (dRef / sunSatVector.getNorm());
        final double rawP   = (pRef * dRatio * dRatio *
                               this.getLightingRatio(s.getPositionVelocity().getPosition(),
                                                     s.getFrame(), s.getDate()));
        /*
         * Spacecraft characteristics effects
         */
        Vector3D u = sunSatVector.normalize();
        Vector3D inSpacecraft = s.getAttitude().getRotation().applyTo(u);
        final double kd = ((1.0 - spacecraft.getAbsorptionCoef(s, inSpacecraft).getNorm()) *
                           (1.0 - spacecraft.getReflectionCoef(s, inSpacecraft).getNorm()));

        final double acceleration =
            (rawP * (1 + kd * 4.0 / 9.0 ) * spacecraft.getRadiationCrossSection(s, inSpacecraft) / s.getMass());
        /*
         * Provide the perturbing acceleration to the derivatives adder
         */
        adder.addXYZAcceleration(acceleration * u.getX(),
                                 acceleration * u.getY(),
                                 acceleration * u.getZ());
    }

    /** 
     * @param position Satellite's position in the selected frame.
     * @param frame Position frame
     * @param date Current date
     * 
     * @return Lighting ratio between zero and one inclusive.
     * 
     * @exception OrekitException Trajectory too low (within
     * equatorial radius)
     */
    public double getLightingRatio(Vector3D position, Frame frame, AbsoluteDate date)
        throws OrekitException
    {
        Vector3D satSunVector =
            sun.getPositionVelocity(date, frame).getPosition().subtract(position);
        /*
         * Earth apparent radius
         */
        final double r = position.getNorm();
        if (r <= equatorialRadius) {
            throw new OrekitException(LOW_TRAJECTORY_MESSAGE, r);
        }

        final double alphaEarth = Math.atan(equatorialRadius / r);
        /*
         * Definition of the Sun's apparent radius
         */
        final double alphaSun = SUN_RADIUS / satSunVector.getNorm();
        /*
         * Retrieve the Sat-Sun / Sat-Central body angle
         */
        final double sunEarthAngle = Vector3D.angle(satSunVector, position.negate());
        /*
         * Is the satellite in complete penumbra ?
         */
        if (sunEarthAngle - alphaEarth + alphaSun < 0.0)
            return 0.0;
        /*
         * Compute a lightning ratio in penumbra
         */
        else if ((sunEarthAngle - alphaEarth + alphaSun >= 0.0) &&
                 (sunEarthAngle - alphaEarth - alphaSun <= 0.0))
        {
            final double alpha1 =
                (sunEarthAngle * sunEarthAngle -
                        (alphaEarth - alphaSun) * (alphaSun + alphaEarth)) / (2 * sunEarthAngle);

            final double alpha2 =
                (sunEarthAngle * sunEarthAngle +
                        (alphaEarth - alphaSun) * (alphaSun + alphaEarth)) / (2 * sunEarthAngle);

            final double P1 = Math.PI * alphaSun * alphaSun -
                alphaSun * alphaSun * Math.acos(alpha1 / alphaSun) +
                alpha1 * Math.sqrt(alphaSun * alphaSun - alpha1 * alpha1);

            final double P2 = alphaEarth * alphaEarth * Math.acos(alpha2 / alphaEarth) -
                alpha2 * Math.sqrt(alphaEarth * alphaEarth - alpha2 * alpha2);

            return ((P1 - P2) / (Math.PI * alphaSun * alphaSun));
        }
        else
            return 1.0;
    }

    /** 
     * @param state current spacecraft state
     * @return Compute Sat-Sun vector in spacecraft state frame.
     * @exception OrekitException if sun position cannot be computed
     */
    private Vector3D getSatSunVector(SpacecraftState state)
        throws OrekitException
    {
        PositionVelocity sunPV = sun.getPositionVelocity(state.getDate(), state.getFrame());
        PositionVelocity satPV = state.getPositionVelocity();
        return sunPV.getPosition().subtract(satPV.getPosition());
    }

    /** {@inheritDoc} */
    public EventDetector[] getEventsDetectors() {
        return new EventDetector[] {
            new UmbraDetector(), new PenumbraDetector()
        };
    }

    /** 
     * This class defines the umbra entry/exit detector.
     */
    private class UmbraDetector
        extends AbstractDetector
    {
        public UmbraDetector() {
            super(60.0, 1.0e-3);
        }


        public int eventOccurred(final SpacecraftState s, final boolean increasing) {
            return RESET_DERIVATIVES;
        }
        /**
         * The G-function is the difference between the
         * Sat-Sun-Sat-Earth angle and the Earth's apparent radius.
         * 
         * @param s  Current state information : date, kinematics, attitude
         * 
         * @return Value of the g function
         * 
         * @exception OrekitException Sun or spacecraft position
         * cannot be computed 
         */
        public double g(SpacecraftState s)
            throws OrekitException
        {
            Vector3D satPos = s.getPositionVelocity().getPosition();
            final double sunEarthAngle = Math.PI - Vector3D.angle(getSatSunVector(s), satPos);
            final double r = satPos.getNorm();
            if (r <= equatorialRadius) {
                throw new OrekitException(LOW_TRAJECTORY_MESSAGE, r);
            }
            else {
                double alphaEarth = equatorialRadius / r;
                return sunEarthAngle - alphaEarth;
            }
        }
    }

    /**
     * This class defines the penumbra entry/exit detector.
     */
    private class PenumbraDetector
        extends AbstractDetector
    {
        public PenumbraDetector() {
            super(60.0, 1.0e-3);
        }


        public int eventOccurred(final SpacecraftState s, final boolean increasing) {
            return RESET_DERIVATIVES;
        }
        /**
         * The G-function is the difference between the Sat-Sun-Sat-Earth angle and
         * the sum of the Earth's and Sun's apparent radius.
         * 
         * @param s Current state information : date, kinematics,
         * attitude
         * 
         * @return Value of the g function
         * 
         * @exception OrekitException Sun or spacecraft position
         * cannot be computed 
         */
        public double g(final SpacecraftState s)
            throws OrekitException
        {
            Vector3D satPos       = s.getPositionVelocity().getPosition();
            Vector3D satSunVector = getSatSunVector(s);
            final double sunEarthAngle  = Math.PI - Vector3D.angle(satSunVector, satPos);
            final double r = satPos.getNorm();
            if (r <= equatorialRadius) {
                throw new OrekitException(LOW_TRAJECTORY_MESSAGE, r);
            }
            else {
                double alphaEarth = equatorialRadius / r;
                double alphaSun   = SUN_RADIUS / satSunVector.getNorm();
                return sunEarthAngle - alphaEarth - alphaSun;
            }
        }
    }

}
